# Justifying Distribution
# Step 03: Make appendix tables and figures
# Last updated: 08/15/2023

# Load packages, functions, and datasets ----------------------------------------------

library(plyr)
library(tidyverse)
library(estimatr)
library(texreg)
library(kableExtra)
library(nnet)
library(lmtest)
library(ggfittext)
library(ggthemes)
library(gridExtra)
library(fixest)


load("cleaned_exp1.RData")
load("cleaned_expR1.RData")
load("cleaned_expR2.RData")
load("cleaned_exp2.RData")

# These two functions make regression tables (texreg) and normal tables (kable)
# into Latex format. Usually, the kable tables need additional formatting in
# Latex to make it into the file seen in the SI.
source("make_texreg.R")
source("make_kable.R")

# This function is used for calculating the AIC
myAIC <- function(data) {
  
  2*(data$k+1) + data$nobs * (log(2*pi* (1-data$r.squared)*data$tss/data$nobs ) + 1)
  
}

options(knitr.table.format = "latex")

# Set path here to save your table and figure output to
# Default is that it creates an output folder in working directory.
dir.create("output")
save_path = "output/"

# *Experiment 1 ------------------------------------------------------------

# Table A1 ----------------------------------------------------

#REFORMATTED IN OVERLEAF AS TWO PAGE TABLE

lma1 <- lm_robust(eval_incumbent ~ paydist_cat + paycity_cat + task, 
                  data = df1, cluster = ResponseId)

tablea1 <- texreg(lma1, 
                  digits = 3, 
                  include.ci = FALSE,
       custom.model.names = c("\\makecell{Incumbent\\\\Evaluations\\\\($-$1 to 1)}"),
       custom.coef.map = list("paydist_cat2" = "District Payoff Categories: $-$\\$25 per capita",
                              "paydist_cat3" = "District Payoff Categories: $-$\\$10 per capita",
                              "paydist_cat4" = "District Payoff Categories: $-$\\$5 per capita",
                              "paydist_cat5" = "District Payoff Categories: $-$\\$1 per capita",
                              "paydist_cat6" = "District Payoff Categories: \\$0 per capita",
                              "paydist_cat7" = "District Payoff Categories: \\$1 per capita",
                              "paydist_cat8" = "District Payoff Categories: \\$5 per capita",
                              "paydist_cat9" = "District Payoff Categories: \\$10 per capita",
                              "paydist_cat10" = "District Payoff Categories: \\$25 per capita",
                              "paydist_cat11" = "District Payoff Categories: \\$50 per capita",
                              "paycity_cat2" = "City Payoff Categories: $-$\\$40 per capita",
                              "paycity_cat3" = "City Payoff Categories: $-$\\$30 per capita",
                              "paycity_cat4" = "City Payoff Categories: $-$\\$25 per capita",
                              "paycity_cat5" = "City Payoff Categories: $-$\\$20 per capita",
                              "paycity_cat6" = "City Payoff Categories: $-$\\$15 per capita",
                              "paycity_cat7" = "City Payoff Categories: $-$\\$10 per capita",
                              "paycity_cat8" = "City Payoff Categories: $-$\\$5 per capita",
                              "paycity_cat9" = "City Payoff Categories: $-$\\$1 per capita",
                              "paycity_cat10" = "City Payoff Categories: $-$\\$0.1 per capita",
                              "paycity_cat11" = "City Payoff Categories: \\$0 per capita",
                              "paycity_cat12" = "City Payoff Categories: \\$0.1 per capita",
                              "paycity_cat13" = "City Payoff Categories: \\$1 per capita",
                              "paycity_cat14" = "City Payoff Categories: \\$5 per capita",
                              "paycity_cat15" = "City Payoff Categories: \\$10 per capita",
                              "paycity_cat16" = "City Payoff Categories: \\$15 per capita",
                              "paycity_cat17" = "City Payoff Categories: \\$20 per capita",
                              "paycity_cat18" = "City Payoff Categories: \\$25 per capita",
                              "paycity_cat19" = "City Payoff Categories: \\$30 per capita",
                              "paycity_cat20" = "City Payoff Categories: \\$40 per capita",
                              "paycity_cat21" = "City Payoff Categories: \\$50 per capita",
                              "task2" = "Vignette 2",
                              "task3" = "Vignette 3",
                              "(Intercept)" = "Constant"),
       caption = "Effect of District and City-Wide Returns on 
       Incumbent Evaluations with Returns as Categorical Variables",
       caption.above = TRUE, 
       include.adjrs = FALSE, 
       include.rmse = FALSE,
       booktabs = T,
       float.pos = "H",
       custom.gof.names = c("R$^2$", "Observations", "Respondents"),
       threeparttable = TRUE,
       custom.note = paste("\\item[\\hspace{-5mm}] %stars.",
                    "\\item[\\hspace{-5mm}] \textit{Note:} Models estimated using ordinary least squares regression, with standard errors 
       clustered by respondent."),
       label = "tab:tablea1")

cat(tablea1, sep = "\n", file = paste0(save_path, "tablea1.tex"))

# Table A2 -----------------------------------------------------

lma2 <- lm_robust(vote100 ~ distatorabovezero +
                    distabovezero +
                    distworsecity + city_pc + dist_pc + task,
                  cluster = ResponseId, data = df1, subset = city > 0)

make_texreg(name = "tablea2",
            model = lma2,
            custom.model.names = c("\\makecell{Vote for Incumbent\\\\vs. Challenger\\\\(0 or 100)}"),
            custom.coef.map = list("distatorabovezero" = "District At Least Breaks Even (District \\geq 0)",
                                   "distabovezero" = "District Benefits (District > 0)",
                                   "dist_pc" = "District Returns Per Capita",
                                   "distworsecity" = "District Worse Off than City",
                                   "city_pc" = "City Returns Per Capita",
                                   "task2" = "Vignette 2",
                                   "task3" = "Vignette 3",
                                   "(Intercept)" = "Constant"),
            custom.note = paste("\\item[\\hspace{-5mm}] %stars.",
                                "\\item[\\hspace{-5mm}] \textit{Note:} Models estimated using ordinary least squares regression, with standard errors 
       clustered by respondent."),
            caption = "Model for Simulated Electoral Tradeoffs when City Per Capita Return is Greater than 0,
             Experiment 1 (Table 2)")

# (Kable) Table 2 -----------------------------------------------------------------
# MANUAL FIXES to table starting from below table, don't refresh.

initial_vote_share_win1 <- 
  lma2$coefficients[1] + lma2$coefficients[2] + lma2$coefficients[3] + 
  lma2$coefficients[5] + lma2$coefficients[6]*4

initial_vote_share_loss1 <- 
  lma2$coefficients[1] + lma2$coefficients[4] +
  lma2$coefficients[5]*1 - lma2$coefficients[6]*1

newlosingmember1 <- 
  lma2$coefficients[1] + lma2$coefficients[2] + lma2$coefficients[4] +
  lma2$coefficients[5]*1

# 6-4, 10-0
dist_pc_levels <- c(3.5, 2.5)

newwinningmember1 <- lma2$coefficients[1] + lma2$coefficients[2] + lma2$coefficients[3] + 
  lma2$coefficients[5]*1 + lma2$coefficients[6]*dist_pc_levels

winners1 <- c(initial_vote_share_win1, newwinningmember1)
losers1 <- c(initial_vote_share_loss1, rep(newlosingmember1, 2))

table2 <- as.data.frame(rbind(round(winners1, 2),
                               round(losers1, 2),
                               round(losers1, 2))) %>%
  `rownames<-`(c("Support in Winning Districts",
                 "Support in Break Even Districts",
                 "Support in Losing Districts"))

## Have to replace the N/A categories to 0.
## I manually turn this into a 7 column table, where I add the district returns 
## for the winners, break even, and losers. I also hand add in the % and the 
## parantheses for the number of districts in each category.
make_kable(table2,
           caption = "Simulated Electoral Tradeoffs from Side Payoffs 
      when City-Wide Per Capita Returns are \\$1, Experiment 1",
           col.names = c("\\makecell{4 Winners\\\\0 Break Even\\\\6 Losers}",
                         "\\makecell{4 Winners\\\\2 Break Even\\\\4 Losers}", 
                         "\\makecell{4 Winners\\\\6 Break Even\\\\0 Losers}"),
           footnote = "\\textit{Note:} This table presents the predicted electoral outcomes estimated from Table A2. 
           We assume a 10-district city passes a policy that has a per capita net return of \\$1 to the city as a whole. 
           By type (winning, break even, losing), all districts are assumed to have the same net return. 
           Districts that come out behind have a per capita net return of $-$\\$1. 
           In column 1, we present the vote share in which 4 districts come out ahead with a per capita net return of \\$4. 
           In column 2, the 4 councilors from the winning districts provide payoffs to 2 losing districts so that they 
           reach a per capita return of \\$0, to create a minimum winning coalition of 6. 
           In column 3, the councilors provide payoffs to 6 losing districts to create a universal coalition of 10.",
           add_header = TRUE)

# (Kable) Table A3 -----------------------------------------------------

lma3 <- df1 %>% 
  select(eval_incumbent, vote, eval_project,
         distatorabovezero, distabovezero,
         cityatorabovezero, cityabovezero,
         distworsecity, city_pc, dist_pc)

tablea3 <- as.data.frame(round(apply(lma3, 2, mean, na.rm = TRUE), digits = 3)) %>% 
  cbind(round(apply(lma3, 2, sd, na.rm = TRUE), digits = 3)) %>% 
  cbind(round(apply(lma3, 2, min, na.rm = TRUE), digits = 3)) %>% 
  cbind(round(apply(lma3, 2, max, na.rm = TRUE), digits = 3)) %>% 
  rownames_to_column() %>% 
  mutate(rowname = case_when(rowname == "eval_incumbent" ~ "Incumbent Evaluation",
         rowname == "eval_project" ~ "Project Evaluation",
         rowname == "vote" ~ "Vote for Incumbent",
         rowname == "distatorabovezero" ~ "District At Least Breaks Even (District \\geq 0)",
         rowname == "distabovezero" ~ "District Benefits (District $>$ 0)",
         rowname == "cityatorabovezero" ~ "City At Least Breaks Even (City \\geq 0)",
         rowname == "cityabovezero" ~ "City Benefits (City $>$ 0)",
         rowname == "distworsecity" ~ "District Worse Off than City",
         rowname == "city_pc" ~ "City Returns Per Capita",
         rowname == "dist_pc" ~ "District Returns Per Capita")) 

make_kable(tablea3, col.names = c("Variable", "Mean", "SD", "Min", "Max"),
           caption = "Summary Statistics for Experiment 1")

# (Kable) Table A4 ------------------------------------------------------------
df1_aic <- df1 %>% 
  mutate(cityabove_125m=city >-12500000) %>% 
  mutate(cityabove_10m=city >-10000000) %>% 
  mutate(cityabove_75m=city >-7500000) %>% 
  mutate(cityabove_625m=city >-6250000) %>% 
  mutate(cityabove_5m=city >-5000000) %>% 
  mutate(cityabove_375m=city >-3750000) %>% 
  mutate(cityabove_25m=city >-2500000) %>% 
  mutate(cityabove_1_25m=city >-1250000) %>% 
  mutate(cityabove_250k=city >-250000) %>% 
  mutate(cityabove_25k=city >-25000) %>% 
  mutate(cityabovezero=city > 0) %>% 
  mutate(cityabove125m=city >12500000) %>% 
  mutate(cityabove10m=city >10000000) %>% 
  mutate(cityabove75m=city >7500000) %>% 
  mutate(cityabove625m=city >6250000) %>% 
  mutate(cityabove5m=city >5000000) %>% 
  mutate(cityabove375m=city >3750000) %>% 
  mutate(cityabove25m=city >2500000) %>% 
  mutate(cityabove1_25m=city >1250000) %>% 
  mutate(cityabove250k=city >250000) %>% 
  mutate(cityabove25k=city >25000) %>% 
  
  mutate(cityatorabove_125m=city >=-12500000) %>% 
  mutate(cityatorabove_10m=city >=-10000000) %>% 
  mutate(cityatorabove_75m=city >=-7500000) %>% 
  mutate(cityatorabove_625m=city >=-6250000) %>% 
  mutate(cityatorabove_5m=city >=-5000000) %>% 
  mutate(cityatorabove_375m=city >=-3750000) %>% 
  mutate(cityatorabove_25m=city >=-2500000) %>% 
  mutate(cityatorabove_1_25m=city >=-1250000) %>% 
  mutate(cityatorabove_250k=city >=-250000) %>% 
  mutate(cityatorabove_25k=city >=-25000) %>% 
  mutate(cityatorabovezero=city >=0) %>% 
  mutate(cityatorabove125m=city >=12500000) %>% 
  mutate(cityatorabove10m=city >=10000000) %>% 
  mutate(cityatorabove75m=city >=7500000) %>% 
  mutate(cityatorabove625m=city >=6250000) %>% 
  mutate(cityatorabove5m=city >=5000000) %>% 
  mutate(cityatorabove375m=city >=3750000) %>% 
  mutate(cityatorabove25m=city >=2500000) %>% 
  mutate(cityatorabove1_25m=city >=1250000) %>% 
  mutate(cityatorabove250k=city >=250000) %>% 
  mutate(cityatorabove25k=city >=25000) %>% 
  
  mutate(distabove_125m=dist > -1250000) %>% 
  mutate(distabove_625k=dist > -625000) %>% 
  mutate(distabove_250k=dist > -250000) %>% 
  mutate(distabove_125k=dist > -125000) %>% 
  mutate(distabove_25k=dist > -25000) %>% 
  mutate(distabovezero=dist > 0) %>% 
  mutate(distabove125m=dist > 1250000) %>% 
  mutate(distabove625k=dist > 625000) %>% 
  mutate(distabove250k=dist > 250000) %>% 
  mutate(distabove125k=dist > 125000) %>% 
  mutate(distabove25k=dist > 25000) %>% 
  
  mutate(distatorabove_125m=dist >= -1250000) %>% 
  mutate(distatorabove_625k=dist >= -625000) %>% 
  mutate(distatorabove_250k=dist >= -250000) %>% 
  mutate(distatorabove_125k=dist >= -125000) %>% 
  mutate(distatorabove_25k=dist >=-25000) %>% 
  mutate(distatorabovezero=dist >=0) %>% 
  mutate(distatorabove125m=dist >= 1250000) %>% 
  mutate(distatorabove625k=dist >= 625000) %>% 
  mutate(distatorabove250k=dist >= 250000) %>% 
  mutate(distatorabove125k=dist >= 125000) %>% 
  mutate(distatorabove25k=dist >= 25000)

cityabovevector <- c(
  "cityabove_10m",
  "cityabove_75m",
  "cityabove_625m",
  "cityabove_5m",
  "cityabove_375m",
  "cityabove_25m",
  "cityabove_1_25m",
  "cityabove_250k",
  "cityabove_25k",
  "cityabovezero", 
  "cityabove25k",
  "cityabove250k",
  "cityabove1_25m",
  "cityabove25m",
  "cityabove375m",
  "cityabove5m",
  "cityabove625m",
  "cityabove75m",
  "cityabove10m")

cityatabovevector <- c(
  "cityatorabove_10m",
  "cityatorabove_75m",
  "cityatorabove_625m",
  "cityatorabove_5m",
  "cityatorabove_375m",
  "cityatorabove_25m",
  "cityatorabove_1_25m",
  "cityatorabove_250k",
  "cityatorabove_25k",
  "cityatorabovezero",
  "cityatorabove25k",
  "cityatorabove250k",
  "cityatorabove1_25m",
  "cityatorabove25m",
  "cityatorabove375m",
  "cityatorabove5m",
  "cityatorabove625m",
  "cityatorabove75m",
  "cityatorabove10m")

distabovevector <- c(
  "distabove_625k",
  "distabove_250k",
  "distabove_125k",
  "distabove_25k",
  "distabovezero",
  "distabove25k",
  "distabove125k",
  "distabove250k",
  "distabove625k")

distatabovevector <- c(
  "distatorabove_625k",
  "distatorabove_250k",
  "distatorabove_125k",
  "distatorabove_25k",
  "distatorabovezero",
  "distatorabove25k",
  "distatorabove125k", 
  "distatorabove250k", 
  "distatorabove625k")

name_rows = c("$-$\\$10M","$-$\\$7.5M","$-$\\$6.25M","$-$\\$5M", "$-$\\$3.75M", "$-$\\$2.5M","$-$\\$1.25M","$-$\\$250K", "$-$\\$25K", 
              "\\$0","\\$25K", "\\$250K","\\$1.25M","\\$2.5M", "\\$3.75M","\\$5M",  "\\$6.25M","\\$7.5M",
              "\\$10M")
name_columns = c("$-$\\$625K",
                 "$-$\\$250K",
                 "$-$\\$125K",
                 "$-$\\$25K",
                 "\\$0",
                 "\\$25K",
                 "\\$125K", 
                 "\\$250K", 
                 "\\$625K")

aic_a4 <-matrix(data=NA,ncol=9, nrow=19)

for(i in 1:length(cityabovevector)){
  for(j in 1:length(distabovevector)){
    model <- lm_robust(
      eval_incumbent ~ !!sym(cityabovevector[i]) + !!sym(cityatabovevector[i]) +
        !!sym(distabovevector[j]) + !!sym(distatabovevector[j])  +
        distworsecity + city_pc + dist_pc, data = df1_aic, cluster = ResponseId)
    aic_a4[i,j] <- myAIC(model)
  }
}

aic_a4

#Get minimum in each row and verify that it's the 0 column for bolding later on
minimum_aic <- apply(aic_a4, 1, FUN = min)
# We notice that for the 13th row, or the city return of $1.25M it's not the 0 district return column.
aic_a4 == minimum_aic

tablea4 <- as.data.frame(aic_a4) %>%
  cbind(name_rows) %>% 
  relocate(name_rows)

make_kable(
  tablea4, 
  col.names = c("", name_columns),
  caption = "Akaike Information Criterion (AIC) by City and District Cutpoints, Experiment 1", 
  footnote = "\\textit{Note:} Cells represent the AIC for each model with the given city and district
cutpoint dummy variables, with the rows representing city cutpoint values and
the columns representing district cutpoint values. The model estimated was 
Incumbent Evaluation = $\\beta_0 + \\beta_1 \\times$District Above Cutpoint$+
  \\beta_2 \\times$District At or Above Cutpoint$+ \\beta_3 \\times$Aggregate 
Above Cutpoint$+ \\beta_4 \\times$Aggregate At or Above Cutpoint$+
  \\beta_5 \\times$District Per Capita Less than City Per Capita$+ \\beta_6 \\times$City Returns Per Capita$+
  \\beta_7 \\times$District Returns Per Capita$+ \\epsilon$. Bold values indicate lowest AIC (i.e. best fitting)
model for each given set of city cutpoints.", column_spec = T, resize = TRUE, refresh = TRUE)


# Table A5 ---------------------------------------------------------------

# Replace \renewcommand{\baselinestretch}{1}% for the top

lma5_1 <- lm_robust(eval_incumbent~ distatorabovezero +
                            distabovezero + cityatorabovezero + cityabovezero +
                            distworsecity + city_pc + dist_pc + 
                            pid7 +hhi + non_white + education + task,
                          cluster = ResponseId, data = df1)
lma5_2 <- lm_robust(eval_incumbent~ distatorabovezero +
                              distabovezero + cityatorabovezero + cityabovezero +
                              distworsecity + city_pc + dist_pc + pid7*distatorabovezero +
                              pid7*distabovezero + task,
                            cluster = ResponseId, data = df1)
lma5_3 <- lm_robust(eval_incumbent~ distatorabovezero +
                              distabovezero + cityatorabovezero + cityabovezero +
                              distworsecity + city_pc + dist_pc + hhi*distatorabovezero +
                              hhi*distabovezero + task,
                            cluster = ResponseId, data = df1)
lma5_4 <- lm_robust(eval_incumbent~ distatorabovezero +
                              distabovezero + cityatorabovezero + cityabovezero +
                              distworsecity + city_pc + dist_pc + non_white*distatorabovezero +
                              non_white*distabovezero + task,
                            cluster = ResponseId, data = df1)
lma5_5 <- lm_robust(eval_incumbent~ distatorabovezero +
                              distabovezero + cityatorabovezero + cityabovezero +
                              distworsecity + city_pc + dist_pc + education*distatorabovezero +
                              education*distabovezero + task,
                            cluster = ResponseId, data = df1)

make_texreg(name = "tablea5",
            model = list(lma5_1, lma5_2, lma5_3, lma5_4, lma5_5),
            custom.model.names = c("\\makecell{Incumbent\\\\Evaluations\\\\($-$1 to 1)}", 
                                   "\\makecell{Incumbent\\\\Evaluations\\\\($-$1 to 1)}", 
                                   "\\makecell{Incumbent\\\\Evaluations\\\\($-$1 to 1)}", 
                                   "\\makecell{Incumbent\\\\Evaluations\\\\($-$1 to 1)}", 
                                   "\\makecell{Incumbent\\\\Evaluations\\\\($-$1 to 1)}"),
            custom.coef.map = list("distatorabovezero" = "District At Least Breaks Even (District \\geq 0)",
                                   "distabovezero" = "District Benefits (District > 0)",
                                   "dist_pc" = "District Returns Per Capita",
                                   "distworsecity" = "District Worse Off than City",
                                   "cityatorabovezero" = "City At Least Breaks Even (City \\geq 0)",
                                   "cityabovezero" = "City Benefits (City > 0)",
                                   "city_pc" = "City Returns Per Capita",
                                   "pid7" = "Party (7-Point Scale)",
                                   "distatorabovezero:pid7" = "District At Least Breaks Even $\\times$ Party",
                                   "distabovezero:pid7" = "District Benefits $\\times$ Party",
                                   "hhi" = "Household Income (24-Point Scale)",
                                   "distatorabovezero:hhi" = "District At Least Breaks Even $\\times$ HHI",
                                   "distabovezero:hhi" = "District Benefits $\\times$ HHI",
                                   "non_white" = "Non-White",
                                   "distatorabovezero:non_white" = "District At Least Breaks Even $\\times$ Non-White",
                                   "distabovezero:non_white" = "District Benefits $\\times$ Non-White",
                                   "education" = "Education (8-Point Scale)",
                                   "distatorabovezero:education" = "District At Least Breaks Even $\\times$ Education",
                                   "distabovezero:education" = "District Benefits $\\times$ Education",
                                   "task2" = "Vignette 2",
                                   "task3" = "Vignette 3",
                                   "(Intercept)" = "Constant"),
            custom.note = paste("\\item[\\hspace{-5mm}] %stars.",
                                "\\item[\\hspace{-5mm}] \\textit{Note:} Dependent variables are listed in each column. 
       Models estimated using ordinary least squares regression, with standard errors clustered by respondent."),
            caption = "Heterogeneous Effect of District and City-Wide Returns on 
       Incumbent Evaluations by Key Demographic Groups, Experiment 1",
            resize_width = TRUE)

# (Kable) Table A6 ---------------------------------------------------------------
full_lrtest <- function(data){
  data_lr <- data %>% #standardize the datasets to conduct a LR test
    filter(!is.na(pid7) & !is.na(hhi) & !is.na(non_white) & !is.na(education)
           & !is.na(region))
  
  multi_mo_city <- multinom(paycity_cat ~ pid7 + hhi + non_white + education + region,
                            data = data_lr,model=TRUE)
  multi_mo_city_basic <- multinom(paycity_cat ~ 1,
                                  data = data_lr,model=TRUE)
  multi_mo_dist <- multinom(paydist_cat ~ pid7 + hhi + non_white + education + region,
                            data = data_lr,model=TRUE)
  multi_mo_dist_basic <- multinom(paydist_cat ~ 1,
                                  data = data_lr,model=TRUE)
  
  lrtest_city <- lmtest::lrtest(multi_mo_city, multi_mo_city_basic)
  lrtest_dist <- lmtest::lrtest(multi_mo_dist, multi_mo_dist_basic)
  
  cityreturn <- c(lrtest_city$Chisq[2], lrtest_city$`Pr(>Chisq)`[2])
  distreturn <- c(lrtest_dist$Chisq[2], lrtest_dist$`Pr(>Chisq)`[2])
  
  output <- as.data.frame(rbind(cityreturn, distreturn)) %>% 
    rowid_to_column() %>% 
    mutate(rowid = case_when(rowid == 1 ~ "City Return Treatment",
                             rowid == 2 ~ "District Return Treatment"))
}

tablea6 <- full_lrtest(df1)

make_kable(tablea6, col.names = c("", "$\\chi^2$",
                    "$p$-value"),
      caption = "Likelihood Ratio Tests, Experiment 1",
  footnote = "\\textit{Note:} We conduct a multinomial logistic regression of the treatment variable
               on covariates party identification, household income, race, education,
               and region of the survey taker to test if our treatments of city and 
               district net return were associated with any key covariates. We then conducted
               a likelihood ratio test of this model with the null model for each treatment variable outcome.
               Each column reflects the likelihood ratio difference between a saturated and null model,
               as well as the $p$-value for rejecting the null hypothesis. Since we fail to reject
               the null for both variables, we are unable to say there is a significant
               association between our treatment variables and these covariates.")

# Table A7 ---------------------------------------------------------------

lma7_1 <- lm_robust(eval_incumbent~ distatorabovezero +
                              distabovezero + cityatorabovezero +
                              distworsecity + city_pc + dist_pc + 
                              distatorabovezero*cityatorabovezero + 
                              distabovezero*cityatorabovezero + 
                              dist_pc*cityatorabovezero + 
                              task,
                            cluster = ResponseId, data = df1)
lma7_2 <- lm_robust(eval_incumbent~ distatorabovezero +
                              distabovezero + cityatorabovezero +
                              distworsecity + city_pc + dist_pc + 
                              distatorabovezero*cityatorabovezero + 
                              distabovezero*cityatorabovezero + 
                              dist_pc*cityatorabovezero + 
                              pid7 + hhi + non_white + education + region+
                              task,
                            cluster = ResponseId, data = df1)

make_texreg(name = "tablea7",
            model = list(lma7_1, lma7_2),
            custom.model.names = c("\\makecell{Incumbent\\\\Evaluations\\\\($-$1 to 1)}", 
                                   "\\makecell{Incumbent\\\\Evaluations\\\\with Covariates\\\\($-$1 to 1)}"),
            custom.coef.map = list("distatorabovezero" = "District At Least Breaks Even (District \\geq 0)",
                                   "distabovezero" = "District Benefits (District > 0)",
                                   "dist_pc" = "District Returns Per Capita",
                                   "distworsecity" = "District Worse Off than City",
                                   "cityatorabovezero" = "City At Least Breaks Even (City \\geq 0)",
                                   "cityabovezero" = "City Benefits (City > 0)",
                                   "city_pc" = "City Returns Per Capita",
                                   "distatorabovezero:cityatorabovezero" = "District At Least Breaks Even $\\times$ City At Least Breaks Even",
                                   "distabovezero:cityatorabovezero" = "District Benefits $\\times$ City At Least Breaks Even",
                                   "cityatorabovezero:dist_pc" = "District Returns Per Capita $\\times$ City At Least Breaks Even",
                                   "task2" = "Vignette 2",
                                   "task3" = "Vignette 3",
                                   "(Intercept)" = "Constant"),
            custom.note = paste("\\item[\\hspace{-5mm}] %stars.",
                                "\\item[\\hspace{-5mm}] \\textit{Note:} Dependent variables are listed in each column. 
       Models estimated using ordinary least squares regression, 
                                with standard errors clustered by respondent. 
                                Covariates included in Column 2 include party
                                identification, household income, 
                                education level, geographic region, and race."), 
            resize_width = TRUE,
            caption = "Effect of District Returns on Incumbent Evaluations by City Outcome, Experiment 1")

# Table A8 ---------------------------------------------------------------

lma8 <- lm_robust(eval_incumbent~ distatorabovezero +
                              distabovezero + cityatorabovezero +
                              distworsecity + city_pc + dist_pc + 
                              big_dist_return +
                              big_dist_return*city_pc + 
                              big_dist_return*dist_pc + 
                              task,
                            cluster = ResponseId, data = df1)

make_texreg(name = "tablea8",
            model = lma8,
            custom.model.names = c("\\makecell{Incumbent\\\\Evaluations\\\\($-$1 to 1)}"),
            custom.coef.map = list("distatorabovezero" = "District At Least Breaks Even (District \\geq 0)",
                                   "distabovezero" = "District Benefits (District > 0)",
                                   "dist_pc" = "District Returns Per Capita",
                                   "distworsecity" = "District Worse Off than City",
                                   "cityatorabovezero" = "City At Least Breaks Even (City \\geq 0)",
                                   "cityabovezero" = "City Benefits (City > 0)",
                                   "city_pc" = "City Returns Per Capita",
                                   "big_dist_return" = "District Return is Substantial and Positive",
                                   "city_pc:big_dist_return" = "Substantial District Return $\\times$ City Returns PC",
                                   "dist_pc:big_dist_return" = "Substantial District Return $\\times$ District Returns PC",
                                   "task2" = "Vignette 2",
                                   "task3" = "Vignette 3",
                                   "(Intercept)" = "Constant"),
            custom.note = paste("\\item[\\hspace{-5mm}] %stars.",
                                "\\item[\\hspace{-5mm}] \\textit{Note:} Dependent variables are listed in each column. 
       Models estimated using ordinary least squares regression, with standard errors clustered by respondent. 
                                Substantial district returns are those that are
                                greater than 10."),
            caption = "Effect of Districts and City Returns when Returns are 
       Substantial and Positive, Experiment 1")


# Table A9 ---------------------------------------------------------------


lma9 <- lm_robust(eval_incumbent~ distatorabovezero +
                              distabovezero + cityatorabovezero + cityabovezero + 
                               city_pc + dist_pc + 
                              task,
                            cluster = ResponseId, data = df1, subset = city_pc > dist_pc)
make_texreg(name = "tablea9",
            model = lma9,
            custom.model.names = c("\\makecell{Incumbent\\\\Evaluations\\\\($-$1 to 1)}"),
            custom.coef.map = list("distatorabovezero" = "District At Least Breaks Even (District \\geq 0)",
                                   "distabovezero" = "District Benefits (District > 0)",
                                   "dist_pc" = "District Returns Per Capita",
                                   "distworsecity" = "District Worse Off than City",
                                   "cityatorabovezero" = "City At Least Breaks Even (City \\geq 0)",
                                   "cityabovezero" = "City Benefits (City > 0)",
                                   "city_pc" = "City Returns Per Capita",
                                   "task2" = "Vignette 2",
                                   "task3" = "Vignette 3",
                                   "(Intercept)" = "Constant"),
            custom.note = paste("\\item[\\hspace{-5mm}] %stars.",
                                "\\item[\\hspace{-5mm}] \\textit{Note:} Dependent variables are listed in each column. 
       Models estimated using ordinary least squares regression, with standard errors clustered by respondent."),
            caption = "Effect of District and City Net Returns when the City Per Capita Return is Greater than the District Per Capita Return, Experiment 1")

# Table A10 ---------------------------------------------------------------
lma10_1 <- lm_robust(eval_incumbent~ distatorabovezero +
                        distabovezero + distworsecity+ cityatorabovezero + cityabovezero +
                        city_pc + dist_pc,
                      cluster = ResponseId, data = df1, subset = task == 1)
lma10_2 <- lm_robust(eval_incumbent~ distatorabovezero + 
                             distabovezero + distworsecity+cityatorabovezero + cityabovezero +
                             city_pc + dist_pc,
                             cluster = ResponseId, data = df1, subset = task == 2)
lma10_3 <- lm_robust(eval_incumbent~ distatorabovezero +
                             distabovezero + distworsecity+ cityatorabovezero + cityabovezero +
                             city_pc + dist_pc,
                             cluster = ResponseId, data = df1, subset = task == 3)

make_texreg(name = "tablea10",
            model = list(lma10_1, lma10_2, lma10_3),
            custom.model.names = c("Round 1", "Round 2", "Round 3"),
            custom.coef.map = list("distatorabovezero" = "District At Least Breaks Even (District \\geq 0)",
                                   "distabovezero" = "District Benefits (District > 0)",
                                   "dist_pc" = "District Returns Per Capita",
                                   "distworsecity" = "District Worse Off than City",
                                   "cityatorabovezero" = "City At Least Breaks Even (City \\geq 0)",
                                   "cityabovezero" = "City Benefits (City > 0)",
                                   "city_pc" = "City Returns Per Capita",
                                   "(Intercept)" = "Constant"),
            custom.note = paste("\\item[\\hspace{-5mm}] %stars.",
                                "\\item[\\hspace{-5mm}] \\textit{Note:} The dependent
                                variable in each column is 
                                Incumbent Evaluation, coded $-$1 to 1. 
                                Models estimated using ordinary
                           least squares regression, with standard errors clustered
                           by respondent."),
            caption = "Effect of District and City-Wide Returns on 
       Incumbent Evaluations by Round, Experiment 1",
            multiple.tasks = "normal")

# *Experiment R1 ----------------------------------------------------------


# Table B1 ----------------------------------------------------------------

table1Rincum <- lm_robust(eval_incumbent~ distatorabovezero +
                            distabovezero + cityatorabovezero + cityabovezero +
                            distworsecity + city_pc + dist_pc + task,
                          cluster = ResponseId, data = df1R1)
table1Rvote <- lm_robust(vote~ distatorabovezero +
                           distabovezero + cityatorabovezero + cityabovezero +
                           distworsecity + city_pc + dist_pc + task,
                         cluster = ResponseId, data = df1R1)
table1Rproj <- lm_robust(eval_project~ distatorabovezero +
                           distabovezero + cityatorabovezero + cityabovezero +
                           distworsecity + city_pc + dist_pc + task,
                         cluster = ResponseId, data = df1R1)

make_texreg(name = "tablea11",
            model = list(table1Rincum, table1Rvote, table1Rproj), 
            custom.model.names = c("\\makecell{Incumbent\\\\Evaluations\\\\($-$1 to 1)}", 
                                   "\\makecell{Vote for Incumbent\\\\vs. Challenger\\\\($-$1 to 1)}",
                                   "\\makecell{Project\\\\Evaluation\\\\($-$1 to 1)}"),
            custom.coef.map = list("distatorabovezero" = "District At Least Breaks Even (District \\geq 0)",
                                   "distabovezero" = "District Benefits (District > 0)",
                                   "dist_pc" = "District Returns Per Capita",
                                   "distworsecity" = "District Worse Off than City",
                                   "cityatorabovezero" = "City At Least Breaks Even (City \\geq 0)",
                                   "cityabovezero" = "City Benefits (City > 0)",
                                   "city_pc" = "City Returns Per Capita",
                                   "task2" = "Vignette 2",
                                   "task3" = "Vignette 3",
                                   "(Intercept)" = "Constant"),
            custom.note = paste("\\item[\\hspace{-5mm}] %stars.",
                                "\\item[\\hspace{-5mm}] \\textit{Note:} Dependent variables are listed in each column. 
       Models estimated using ordinary least squares regression, with standard errors clustered by respondent."),
            caption = "Effect of District and City-Wide Returns on 
        Evaluations, Experiment R1", resize_width = TRUE)


# Table B2 ----------------------------------------------------------------

lma12 <- lm_robust(vote100 ~ distatorabovezero +
                    distabovezero +
                    distworsecity + city_pc + dist_pc + task,
                  cluster = ResponseId, data = df1R1, subset = city > 0)

make_texreg(name = "tablea12",
            model = lma12,
            custom.model.names = c("\\makecell{Vote for Incumbent\\\\vs. Challenger\\\\(0 or 100)}"),
            custom.coef.map = list("distatorabovezero" = "District At Least Breaks Even (District \\geq 0)",
                                   "distabovezero" = "District Benefits (District > 0)",
                                   "dist_pc" = "District Returns Per Capita",
                                   "distworsecity" = "District Worse Off than City",
                                   "city_pc" = "City Returns Per Capita",
                                   "task2" = "Vignette 2",
                                   "task3" = "Vignette 3",
                                   "(Intercept)" = "Constant"),
            custom.note = paste("\\item[\\hspace{-5mm}] %stars.",
                                "\\item[\\hspace{-5mm}] \\textit{Note:} Models estimated using ordinary least squares regression, 
                           with standard errors clustered by respondent."),
            caption = "Model for Simulated Electoral Tradeoffs when City Per Capita Return is Greater than 0,
             Experiment R1 (Table B3)")

# (Kable) Table B3 ----------------------------------------------------------------
# MANUAL - see similar process to Table 2 fixes above.

initial_vote_share_win2 <- 
  lma12$coefficients[1] + lma12$coefficients[2] + lma12$coefficients[3] + 
  lma12$coefficients[5] + lma12$coefficients[6]*4

initial_vote_share_loss2 <- 
  lma12$coefficients[1] + lma12$coefficients[4] +
  lma12$coefficients[5]*1 - lma12$coefficients[6]*1

newlosingmember2 <- 
  lma12$coefficients[1] + lma12$coefficients[2] + lma12$coefficients[4] +
  lma12$coefficients[5]*1

dist_pc_levels <- c(3.5, 2.5)

newwinningmember2 <- lma12$coefficients[1] + lma12$coefficients[2] + lma12$coefficients[3] + 
  lma12$coefficients[5]*1 + lma12$coefficients[6]*dist_pc_levels

winners2 <- c(initial_vote_share_win2, newwinningmember2)
losers2 <- c(initial_vote_share_loss2, rep(newlosingmember2, 2))

tablea13 <- as.data.frame(rbind(round(winners2, 2),
                              round(losers2, 2),
                              round(losers2, 2))) %>%
  `rownames<-`(c("Support in Winning Districts",
                 "Support in Break Even Districts",
                 "Support in Losing Districts"))

make_kable(tablea13,
           caption = "Simulated Electoral Tradeoffs from Side Payoffs 
      when City-Wide Per Capita Returns are \\$1, Experiment R1",
           col.names = c("\\makecell{4 Winners\\\\0 Break Even\\\\6 Losers}",
                         "\\makecell{4 Winners\\\\2 Break Even\\\\4 Losers}", 
                         "\\makecell{4 Winners\\\\6 Break Even\\\\0 Losers}"),
           footnote = "\\textit{Note:} This table presents the predicted electoral outcomes estimated from Table B2. 
           We assume a 10-district city passes a policy that has a per capita net return of \\$1 to the city as a whole. 
           By type (winning, break even, losing), all districts are assumed to have the same net return. 
           Districts that come out behind have a per capita net return of $-$\\$1. 
           In column 1, we present the vote share in which 4 districts come out ahead with a per capita net return of \\$4. 
           In column 2, the 4 councilors from the winning districts provide payoffs to 2 losing districts so that they 
           reach a per capita return of \\$0, to create a minimum winning coalition of 6. 
           In column 3, the councilors provide payoffs to 6 losing districts to create a universal coalition of 10.",
           add_header = TRUE)

# (Kable) Table B4 ----------------------------------------------------------------

df1R1_aic <- df1R1 %>% 
  mutate(cityabove_125m=city >-12500000) %>% 
  mutate(cityabove_10m=city >-10000000) %>% 
  mutate(cityabove_75m=city >-7500000) %>% 
  mutate(cityabove_625m=city >-6250000) %>% 
  mutate(cityabove_5m=city >-5000000) %>% 
  mutate(cityabove_375m=city >-3750000) %>% 
  mutate(cityabove_25m=city >-2500000) %>% 
  mutate(cityabove_1_25m=city >-1250000) %>% 
  mutate(cityabove_250k=city >-250000) %>% 
  mutate(cityabove_25k=city >-25000) %>% 
  mutate(cityabovezero=city > 0) %>% 
  mutate(cityabove125m=city >12500000) %>% 
  mutate(cityabove10m=city >10000000) %>% 
  mutate(cityabove75m=city >7500000) %>% 
  mutate(cityabove625m=city >6250000) %>% 
  mutate(cityabove5m=city >5000000) %>% 
  mutate(cityabove375m=city >3750000) %>% 
  mutate(cityabove25m=city >2500000) %>% 
  mutate(cityabove1_25m=city >1250000) %>% 
  mutate(cityabove250k=city >250000) %>% 
  mutate(cityabove25k=city >25000) %>% 
  
  mutate(cityatorabove_125m=city >=-12500000) %>% 
  mutate(cityatorabove_10m=city >=-10000000) %>% 
  mutate(cityatorabove_75m=city >=-7500000) %>% 
  mutate(cityatorabove_625m=city >=-6250000) %>% 
  mutate(cityatorabove_5m=city >=-5000000) %>% 
  mutate(cityatorabove_375m=city >=-3750000) %>% 
  mutate(cityatorabove_25m=city >=-2500000) %>% 
  mutate(cityatorabove_1_25m=city >=-1250000) %>% 
  mutate(cityatorabove_250k=city >=-250000) %>% 
  mutate(cityatorabove_25k=city >=-25000) %>% 
  mutate(cityatorabovezero=city >=0) %>% 
  mutate(cityatorabove125m=city >=12500000) %>% 
  mutate(cityatorabove10m=city >=10000000) %>% 
  mutate(cityatorabove75m=city >=7500000) %>% 
  mutate(cityatorabove625m=city >=6250000) %>% 
  mutate(cityatorabove5m=city >=5000000) %>% 
  mutate(cityatorabove375m=city >=3750000) %>% 
  mutate(cityatorabove25m=city >=2500000) %>% 
  mutate(cityatorabove1_25m=city >=1250000) %>% 
  mutate(cityatorabove250k=city >=250000) %>% 
  mutate(cityatorabove25k=city >=25000) %>% 
  
  mutate(distabove_125m=dist > -1250000) %>% 
  mutate(distabove_625k=dist > -625000) %>% 
  mutate(distabove_250k=dist > -250000) %>% 
  mutate(distabove_125k=dist > -125000) %>% 
  mutate(distabove_25k=dist > -25000) %>% 
  mutate(distabovezero=dist > 0) %>% 
  mutate(distabove125m=dist > 1250000) %>% 
  mutate(distabove625k=dist > 625000) %>% 
  mutate(distabove250k=dist > 250000) %>% 
  mutate(distabove125k=dist > 125000) %>% 
  mutate(distabove25k=dist > 25000) %>% 
  
  mutate(distatorabove_125m=dist >= -1250000) %>% 
  mutate(distatorabove_625k=dist >= -625000) %>% 
  mutate(distatorabove_250k=dist >= -250000) %>% 
  mutate(distatorabove_125k=dist >= -125000) %>% 
  mutate(distatorabove_25k=dist >=-25000) %>% 
  mutate(distatorabovezero=dist >=0) %>% 
  mutate(distatorabove125m=dist >= 1250000) %>% 
  mutate(distatorabove625k=dist >= 625000) %>% 
  mutate(distatorabove250k=dist >= 250000) %>% 
  mutate(distatorabove125k=dist >= 125000) %>% 
  mutate(distatorabove25k=dist >= 25000)

aic_a14 <-matrix(data=NA,ncol=9, nrow=19)

for(i in 1:length(cityabovevector)){
  for(j in 1:length(distabovevector)){
    model <- lm_robust(
      eval_incumbent ~ !!sym(cityabovevector[i]) + !!sym(cityatabovevector[i]) +
        !!sym(distabovevector[j]) + !!sym(distatabovevector[j])  +
        distworsecity + city_pc + dist_pc, data = df1R1_aic, cluster = ResponseId)
    aic_a14[i,j] <- myAIC(model)
  }
}

aic_a14

#Get minimum in each row and verify that it's the 0 column
minimum_aic <- apply(aic_a14, 1, FUN = min)
aic_a14 == minimum_aic

tablea14 <- as.data.frame(aic_a14) %>%
  cbind(name_rows) %>% 
  relocate(name_rows)

make_kable(tablea14,
           caption = "Akaike Information Criterion (AIC) by City and District Cutpoints, Experiment R1",
           col.names = c("", name_columns), resize = TRUE,
           footnote = "\\textit{Note:} Cells represent the AIC for each model with the given city and district
               cutpoint dummy variables, with the rows representing city cutpoint values and
               the columns representing district cutpoint values. The model estimated was 
Incumbent Evaluation = $\\beta_0 + \\beta_1 \\times$District Above Cutpoint$+
  \\beta_2 \\times$District At or Above Cutpoint$+ \\beta_3 \\times$Aggregate 
Above Cutpoint$+ \\beta_4 \\times$Aggregate At or Above Cutpoint$+
  \\beta_5 \\times$District Per Capita Less than City Per Capita$+ \\beta_6 \\times$City Returns Per Capita$+
  \\beta_7 \\times$District Returns Per Capita$+ \\epsilon$. Bold values indicate lowest AIC (i.e. best fitting)
               model for each given set of city cutpoints.",
           column_spec = T)


# (Kable) Table B5 -----------------------------------------------------


tablea15 <- full_lrtest(df1R1)
make_kable(tablea15, col.names = c("", "$\\chi^2$",
                                   "$p$-value"),
           caption = "Likelihood Ratio Tests, Experiment R1",
           footnote = "\\textit{Note:} We conduct a multinomial logistic regression of the treatment variable
               on covariates party identification, household income, race, education,
               and region of the survey taker to test if our treatments of city and 
               district net return were associated with any key covariates. We then conducted
               a likelihood ratio test of this model with the null model for each treatment variable outcome.
               Each column reflects the likelihood ratio difference between a saturated and null model,
               as well as the $p$-value for rejecting the null hypothesis. Since we fail to reject
               the null for both variables, we are unable to say there is a significant
               association between our treatment variables and these covariates.")
# *Experiment R2 -----------------------------------------------------------


# Table C1 ----------------------------------------------------------------
table1R2incum <- lm_robust(eval_incumbent~ distatorabovezero +
                             distabovezero +
                             distworsecity + city_pc + dist_pc + task,
                           cluster = ResponseId, data = df1R2, subset = aggregate_order == "County")
table1R2vote <- lm_robust(vote~ distatorabovezero +
                            distabovezero +
                            distworsecity + city_pc + dist_pc + task,
                          cluster = ResponseId, data = df1R2, subset = aggregate_order == "County")
table1R2proj <- lm_robust(eval_project~ distatorabovezero +
                            distabovezero  +
                            distworsecity + city_pc + dist_pc + task,
                          cluster = ResponseId, data = df1R2, subset = aggregate_order == "County")

make_texreg(name = "tablea16",
            model = list(table1R2incum, table1R2vote, table1R2proj),
            custom.model.names = c("\\makecell{Incumbent\\\\Evaluations\\\\(-1 to 1)}", 
                                   "\\makecell{Vote for Incumbent\\\\vs. Challenger\\\\(-1 to 1)}",
                                   "\\makecell{Project\\\\Evaluation\\\\(-1 to 1)}"),
            custom.coef.map = list("distatorabovezero" = "District At Least Breaks Even (District $\\geq 0$)",
                                   "distabovezero" = "District Benefits (District > 0)",
                                   "distworsecity" = "District Worse Off than County",
                                   "dist_pc" = "District Returns Per Capita",
                                   "city_pc" = "County Returns Per Capita",
                                   "task2" = "Vignette 2",
                                   "task3" = "Vignette 3",
                                   "task4" = "Vignette 4",
                                   "task5" = "Vignette 5",
                                   "(Intercept)" = "Constant"),
            custom.note = paste("\\item[\\hspace{-5mm}] %stars.",
                                "\\item[\\hspace{-5mm}] \\textit{Note:} Dependent variables are listed in each column. 
           Models estimated using ordinary least squares regression, with standard errors clustered by respondent."),
                caption = "Effect of District and County-Wide Returns on 
            Evaluations, Experiment R2", resize_width = TRUE)


# Table C2 ----------------------------------------------------------------


table1R2incum <- lm_robust(eval_incumbent~ distatorabovezero +
                             distabovezero +
                             distworsecity + city_pc + dist_pc + task,
                           cluster = ResponseId, data = df1R2, subset = aggregate_order == "State")
table1R2vote <- lm_robust(vote~ distatorabovezero +
                            distabovezero +
                            distworsecity + city_pc + dist_pc + task,
                          cluster = ResponseId, data = df1R2, subset = aggregate_order == "State")
table1R2proj <- lm_robust(eval_project~ distatorabovezero +
                            distabovezero  +
                            distworsecity + city_pc + dist_pc + task,
                          cluster = ResponseId, data = df1R2, subset = aggregate_order == "State")


make_texreg(name = "tablea17",
            model = list(table1R2incum, table1R2vote, table1R2proj),
            custom.model.names = c("\\makecell{Incumbent\\\\Evaluations\\\\(-1 to 1)}", 
                                   "\\makecell{Vote for Incumbent\\\\vs. Challenger\\\\(-1 to 1)}",
                                   "\\makecell{Project\\\\Evaluation\\\\(-1 to 1)}"),
            custom.coef.map = list("distatorabovezero" = "District At Least Breaks Even (District $\\geq 0$)",
                                   "distabovezero" = "District Benefits (District > 0)",
                                   "distworsecity" = "District Worse Off than State",
                                   "dist_pc" = "District Returns Per Capita",
                                   "city_pc" = "State Returns Per Capita",
                                   "task2" = "Vignette 2",
                                   "task3" = "Vignette 3",
                                   "task4" = "Vignette 4",
                                   "task5" = "Vignette 5",
                                   "(Intercept)" = "Constant"),
            custom.note = paste("\\item[\\hspace{-5mm}] %stars.",
                                "\\item[\\hspace{-5mm}] \\textit{Note:} Dependent variables are listed in each column. 
       Models estimated using ordinary least squares regression, with standard errors clustered by respondent."),
            caption = "Effect of District and State-Wide Returns on 
            Evaluations, Experiment R2", resize_width = TRUE)


# Table C3 ----------------------------------------------------------------

# \renewcommand{\baselinestretch}{1}%

table1R2incum <- lm_robust(eval_incumbent~distatorabovezero*aggregate_order +
                             distabovezero*aggregate_order +
                             distworsecity*aggregate_order + city_pc*aggregate_order + dist_pc*aggregate_order + task*aggregate_order,
                           cluster = ResponseId, data = df1R2)

table1R2vote <- lm_robust(vote~ distatorabovezero*aggregate_order +
                            distabovezero*aggregate_order +
                            distworsecity*aggregate_order + city_pc*aggregate_order + dist_pc*aggregate_order + task*aggregate_order,
                          cluster = ResponseId, data = df1R2)
table1R2proj <- lm_robust(eval_project~ distatorabovezero*aggregate_order +
                            distabovezero*aggregate_order +
                            distworsecity*aggregate_order + city_pc*aggregate_order + dist_pc*aggregate_order + task*aggregate_order,
                          cluster = ResponseId, data = df1R2)

make_texreg(name = "tablea18",
            model = list(table1R2incum, table1R2vote, table1R2proj),
            custom.model.names = c("\\makecell{Incumbent\\\\Evaluations\\\\(-1 to 1)}", 
                                   "\\makecell{Vote for Incumbent\\\\vs. Challenger\\\\(-1 to 1)}",
                                   "\\makecell{Project\\\\Evaluation\\\\(-1 to 1)}"),
            custom.coef.map = list("distatorabovezero" = "District At Least Breaks Even (District $\\geq 0$)",
                                   "distabovezero" = "District Benefits (District > 0)",
                                   "distworsecity" = "District Worse Off than Aggregate",
                                   "dist_pc" = "District Returns Per Capita",
                                   "city_pc" = "Aggregate Returns Per Capita",
                                   "task2" = "Vignette 2",
                                   "task3" = "Vignette 3",
                                   "task4" = "Vignette 4",
                                   "task5" = "Vignette 5",
                                   "aggregate_orderState" = "State",
                                   "distatorabovezero:aggregate_orderState" = "State x District At Least Breaks Even",
                                   "aggregate_orderState:distabovezero" = "State x District Benefits",
                                   "aggregate_orderState:distworsecity" = "State x District Worse Off than Aggregate",
                                   "aggregate_orderState:dist_pc" = "State x District Returns Per Capita",
                                   "aggregate_orderState:city_pc" = "State x Aggregate Returns Per Capita",
                                   "aggregate_orderState:task2" = "State x Vignette 2",
                                   "aggregate_orderState:task3"= "State x Vignette 3",
                                   "aggregate_orderState:task4" = "State x Vignette 4",
                                   "aggregate_orderState:task5" = "State x Vignette 5",
                                   "(Intercept)" = "Constant"),
            custom.note = paste("\\item[\\hspace{-5mm}] %stars.",
                                "\\item[\\hspace{-5mm}] \\textit{Note:} Dependent variables are listed in each column. 
       Models estimated using ordinary least squares regression, with standard errors clustered by respondent."),
            caption = "Effect of District and Aggregate-Wide Returns on 
            Evaluations with State Interactions, Experiment R2",
            refresh = TRUE, scalebox = 0.85)


# Table C4 ----------------------------------------------------------------


table1R2incum <- feols(eval_incumbent~ distatorabovezero +
                         distabovezero +
                         distworsecity + city_pc + dist_pc + task | ResponseId, df1R2, 
                       vcov = ~ResponseId)

table1R2vote <- feols(vote~ distatorabovezero +
                        distabovezero +
                        distworsecity + city_pc + dist_pc + task | ResponseId, df1R2, 
                      vcov = ~ResponseId)

table1R2proj <- feols(eval_project~ distatorabovezero +
                        distabovezero +
                        distworsecity + city_pc + dist_pc + task | ResponseId, df1R2, 
                      vcov = ~ResponseId)


make_texreg(name = "tablea19",
            model = list(table1R2incum, table1R2vote, table1R2proj),
            multiple.tasks = "fe",
            custom.model.names = c("\\makecell{Incumbent\\\\Evaluations\\\\(-1 to 1)}", 
                                   "\\makecell{Vote for Incumbent\\\\vs. Challenger\\\\(-1 to 1)}",
                                   "\\makecell{Project\\\\Evaluation\\\\(-1 to 1)}"),
            custom.coef.map = list("distatorabovezero" = "District At Least Breaks Even (District $\\geq 0$)",
                                   "distabovezero" = "District Benefits (District > 0)",
                                   "distworsecity" = "District Worse Off than Aggregate",
                                   "dist_pc" = "District Returns Per Capita",
                                   "city_pc" = "Aggregate Returns Per Capita",
                                   "task2" = "Vignette 2",
                                   "task3" = "Vignette 3",
                                   "task4" = "Vignette 4",
                                   "task5" = "Vignette 5",
                                   "(Intercept)" = "Constant"),
            custom.note = paste("\\item[\\hspace{-5mm}] %stars.",
                                "\\item[\\hspace{-5mm}] \\textit{Note:} Dependent variables are listed in each column. 
       Models estimated using ordinary least squares regression, with standard errors clustered by respondent."),
            caption = "Effect of District and Aggregate-Wide Returns on Evaluations with Fixed Effects by Individual as Unit (Within Person Analysis), Experiment R2",
            resize_width = TRUE)


# Table C5 ----------------------------------------------------------------


# 3. Only Task 1 or 2

table1R2incum <- lm_robust(eval_incumbent~ distatorabovezero +
                             distabovezero +
                             distworsecity + city_pc + dist_pc + task,
                           cluster = ResponseId, data = df1R2, subset = task == 1|task==2)
table1R2vote <- lm_robust(vote~ distatorabovezero +
                            distabovezero +
                            distworsecity + city_pc + dist_pc + task,
                          cluster = ResponseId, data = df1R2, subset = task == 1|task==2)
table1R2proj <- lm_robust(eval_project~ distatorabovezero +
                            distabovezero  +
                            distworsecity + city_pc + dist_pc + task,
                          cluster = ResponseId, data = df1R2, subset = task == 1|task==2)

make_texreg(name = "tablea20",
            model = list(table1R2incum, table1R2vote, table1R2proj),
            custom.model.names = c("\\makecell{Incumbent\\\\Evaluations\\\\(-1 to 1)}", 
                                   "\\makecell{Vote for Incumbent\\\\vs. Challenger\\\\(-1 to 1)}",
                                   "\\makecell{Project\\\\Evaluation\\\\(-1 to 1)}"),
            custom.coef.map = list("distatorabovezero" = "District At Least Breaks Even (District $\\geq 0$)",
                                   "distabovezero" = "District Benefits (District > 0)",
                                   "distworsecity" = "District Worse Off than Aggregate",
                                   "dist_pc" = "District Returns Per Capita",
                                   "city_pc" = "Aggregate Returns Per Capita",
                                   "task2" = "Vignette 2",
                                   "task3" = "Vignette 3",
                                   "task4" = "Vignette 4",
                                   "task5" = "Vignette 5",
                                   "(Intercept)" = "Constant"),
            custom.note = paste("\\item[\\hspace{-5mm}] %stars.",
                                "\\item[\\hspace{-5mm}] \\textit{Note:} Dependent variables are listed in each column. 
       Models estimated using ordinary least squares regression, with standard errors clustered by respondent."),
            caption = "Effect of District and Aggregate-Wide Returns on Evaluations Subsetted for Vignettes 1 and 2, Experiment R2")


# (Kable) Table C6 ----------------------------------------------------------------
dfR2_aic <- df1R2 %>% 
  mutate(cityabove_5=city_pc  > 5) %>% 
  mutate(cityabove_10=city_pc > 10) %>% 
  mutate(cityabove_15=city_pc > 15) %>% 
  mutate(cityabove_25=city_pc > 25) %>% 
  mutate(cityabove_50=city_pc > 50) %>% 
  
  
  mutate(cityatabove_5=city_pc  >= 5) %>% 
  mutate(cityatabove_10=city_pc >= 10) %>% 
  mutate(cityatabove_15=city_pc >= 15) %>% 
  mutate(cityatabove_25=city_pc >= 25) %>% 
  mutate(cityatabove_50=city_pc >= 50) %>% 
  
  mutate(distabove_55=dist_pc > -55) %>% 
  mutate(distabove_50=dist_pc > -50) %>% 
  mutate(distabove_45=dist_pc > -45) %>% 
  mutate(distabove_30=dist_pc > -30) %>% 
  mutate(distabove_25=dist_pc > -25) %>% 
  mutate(distabove_20=dist_pc > -20) %>% 
  mutate(distabove_15=dist_pc > -15) %>% 
  mutate(distabove_10=dist_pc > -10) %>% 
  mutate(distabove_5=dist_pc > -5) %>% 
  mutate(distabove_2=dist_pc > -2) %>% 
  mutate(distabove_0.5=dist_pc > -0.5) %>% 
  mutate(distabove_0.1=dist_pc > -0.1) %>% 
  mutate(distabovezero=dist_pc > 0) %>%   
  mutate(distabove55=dist_pc > 55) %>% 
  mutate(distabove50=dist_pc > 50) %>% 
  mutate(distabove45=dist_pc > 45) %>% 
  mutate(distabove30=dist_pc > 30) %>% 
  mutate(distabove25=dist_pc > 25) %>% 
  mutate(distabove20=dist_pc > 20) %>% 
  mutate(distabove15=dist_pc > 15) %>% 
  mutate(distabove10=dist_pc > 10) %>% 
  mutate(distabove5=dist_pc > 5) %>% 
  mutate(distabove2=dist_pc > 2) %>% 
  mutate(distabove0.5=dist_pc > 0.5) %>% 
  mutate(distabove0.1=dist_pc > 0.1) %>% 
  
  mutate(distatabove_55=dist_pc >= -55) %>% 
  mutate(distatabove_50=dist_pc >= -50) %>% 
  mutate(distatabove_45=dist_pc >= -45) %>% 
  mutate(distatabove_30=dist_pc >= -30) %>% 
  mutate(distatabove_25=dist_pc >= -25) %>% 
  mutate(distatabove_20=dist_pc >= -20) %>% 
  mutate(distatabove_15=dist_pc >= -15) %>% 
  mutate(distatabove_10=dist_pc >= -10) %>% 
  mutate(distatabove_5=dist_pc >= -5) %>% 
  mutate(distatabove_2=dist_pc >= -2) %>% 
  mutate(distatabove_0.5=dist_pc >= -0.5) %>% 
  mutate(distatabove_0.1=dist_pc >= -0.1) %>% 
  mutate(distatabovezero=dist_pc >= 0) %>%   
  mutate(distatabove55=dist_pc >= 55) %>% 
  mutate(distatabove50=dist_pc >= 50) %>% 
  mutate(distatabove45=dist_pc >= 45) %>% 
  mutate(distatabove30=dist_pc >= 30) %>% 
  mutate(distatabove25=dist_pc >= 25) %>% 
  mutate(distatabove20=dist_pc >= 20) %>% 
  mutate(distatabove15=dist_pc >= 15) %>% 
  mutate(distatabove10=dist_pc >= 10) %>% 
  mutate(distatabove5=dist_pc >= 5) %>% 
  mutate(distatabove2=dist_pc >= 2) %>% 
  mutate(distatabove0.5=dist_pc >= 0.5) %>% 
  mutate(distatabove0.1=dist_pc >= 0.1)


cityabovevector <- c(
  "cityabove_5",
  "cityabove_10",
  "cityabove_15",
  "cityabove_25",
  "cityabove_50")


cityatabovevector <- c(
  "cityatabove_5",
  "cityatabove_10",
  "cityatabove_15",
  "cityatabove_25",
  "cityatabove_50")

distabovevector <- c(
  "distabove_55",
  "distabove_50",
  "distabove_45",
  "distabove_30",
  "distabove_25",
  "distabove_20",
  "distabove_15",
  "distabove_10",
  "distabove_5",
  "distabove_2",
  "distabove_0.5",
  "distabove_0.1",
  "distabovezero",
  "distabove0.1",
  "distabove0.5",
  "distabove2",
  "distabove5",
  "distabove10",
  "distabove15",
  "distabove20",
  "distabove25",
  "distabove30",
  "distabove45",
  "distabove50",
  "distabove55")

distatabovevector <- c(
  "distatabove_55",
  "distatabove_50",
  "distatabove_45",
  "distatabove_30",
  "distatabove_25",
  "distatabove_20",
  "distatabove_15",
  "distatabove_10",
  "distatabove_5",
  "distatabove_2",
  "distatabove_0.5",
  "distatabove_0.1",
  "distatabovezero",
  "distatabove0.1",
  "distatabove0.5",
  "distatabove2",
  "distatabove5",
  "distatabove10",
  "distatabove15",
  "distatabove20",
  "distatabove25",
  "distatabove30",
  "distatabove45",
  "distatabove50",
  "distatabove55")

name_rows = c("$-$\\$55","$-$\\$50","$-$\\$45","$-$\\$30", "$-$\\$25", "$-$\\$20","$-$\\$15","$-$\\$10", "$-$\\$5", "$-$\\$2","$-$\\$0.5","$-$\\$0.1",
              "\\$0","\\$0.1", "\\$0.5","\\$2","\\$5", "\\$10","\\$15",  "\\$20","\\$25",
              "\\$30", "\\$45", "\\$50","\\$55")

name_columns = c("\\$5",
                 "\\$10",
                 "\\$15",
                 "\\$25",
                 "\\$50")



aic_c6 <-matrix(data=NA,ncol=5, nrow=25)

for(i in 1:length(cityabovevector)){
  for(j in 1:length(distabovevector)){
    model <- lm_robust(
      eval_incumbent ~ !!sym(cityabovevector[i]) + !!sym(cityatabovevector[i]) +
        !!sym(distabovevector[j]) + !!sym(distatabovevector[j])  +
        distworsecity + city_pc+ dist_pc, data = dfR2_aic, cluster = ResponseId)
    aic_c6[j,i] <- myAIC(model)
  }
}

aic_c6

test_aic_c6 <- t(aic_c6)

#Get minimum in each row and verify that it's the 0 column
minimum_aic <- apply(test_aic_c6, 1, FUN = min)
test_aic_c6 == minimum_aic

tablea21 <- as.data.frame(aic_c6) %>%
  cbind(name_rows) %>% 
  relocate(name_rows)

make_kable(
  tablea21, 
  col.names = c("", name_columns),
  caption = "Akaike Information Criterion (AIC) by City and District Cutpoints, Experiment R2", 
  footnote = "\\textit{Note:} Cells represent the AIC for each model with the given aggregate and district per capita 
cutpoint dummy variables. Unlike previous models, the \\textit{columns} represent the aggregate cutpoint values and
the \\textit{rows} represent the district cutpoint values. The model estimated was 
Incumbent Evaluation = $\\beta_0 + \\beta_1 \\times$District Above Cutpoint$+
  \\beta_2 \\times$District At or Above Cutpoint$+ \\beta_3 \\times$Aggregate 
Above Cutpoint$+ \\beta_4 \\times$Aggregate At or Above Cutpoint$+
  \\beta_5 \\times$District Per Capita Less than City Per Capita$+ \\beta_6 \\times$City Returns Per Capita$+
  \\beta_7 \\times$District Returns Per Capita$+ \\epsilon$. Bold values indicate lowest AIC (i.e. best fitting)
model for each given set of city cutpoints.", column_spec = FALSE)

#\$0 & \textbf{16910.04} & \textbf{16909.87} & \textbf{16908.79} & \textbf{16908.54} & \textbf{16908.72}\\

# (Kable) Table C7 --------------------------------------------------------

tablec7 <- full_lrtest(df1R2)

make_kable(tablec7, col.names = c("", "$\\chi^2$",
                                  "$p$-value"),
           caption = "Likelihood Ratio Tests, Experiment R2",
           footnote = "\\textit{Note:} We conduct a multinomial logistic regression of the treatment variable
               on covariates party identification, household income, race, education,
               and region of the survey taker to test if our treatments of city and 
               district net return were associated with any key covariates. We then conducted
               a likelihood ratio test of this model with the null model for each treatment variable outcome.
               Each column reflects the likelihood ratio difference between a saturated and null model,
               as well as the $p$-value for rejecting the null hypothesis. Since we fail to reject
               the null for both variables, we are unable to say there is a significant
               association between our treatment variables and these covariates.")


# *Experiment 2 -----------------------------------------------------------
# (Kable) Table D1 ----------------------------------------------------------------


lma22 <- df2 %>% 
  select(councilor_thermometer, vote_num, 
         approve_num, project_thermometer, 
         distatorabovezero, distabovezero,
         dist_worse_than_city, city_pc, dist_pc,
         T_generic, T_baddealdist_notgermane,
         T_baddealdist_germane, T_notfairshare_notgermane,
         T_notfairshare_germane)

tablea22 <- as.data.frame(round(apply(lma22, 2, mean, na.rm = TRUE), digits = 3)) %>% 
  cbind(round(apply(lma22, 2, sd, na.rm = TRUE), digits = 3)) %>% 
  cbind(round(apply(lma22, 2, min, na.rm = TRUE), digits = 3)) %>% 
  cbind(round(apply(lma22, 2, max, na.rm = TRUE), digits = 3)) %>% 
  rownames_to_column() %>% 
  mutate(rowname = case_when(rowname == "councilor_thermometer" ~ "Incumbent Evaluation",
                             rowname == "vote_num" ~ "Vote for Incumbent",
                             rowname == "approve_num" ~ "Approval of Project",
                             rowname == "project_thermometer" ~ "Project Evaluation",
                             rowname == "distatorabovezero" ~ "District At Least Breaks Even (District \\geq 0)",
                             rowname == "distabovezero" ~ "District Benefits (District $>$ 0)",
                             rowname == "cityatorabovezero" ~ "City At Least Breaks Even (City \\geq 0)",
                             rowname == "cityabovezero" ~ "City Benefits (City $>$ 0)",
                             rowname == "dist_worse_than_city" ~ "District Worse Off than City",
                             rowname == "city_pc" ~ "City Returns Per Capita",
                             rowname == "dist_pc" ~ "District Returns Per Capita",
                             rowname == "T_generic" ~ "Generic",
                             rowname == "T_baddealdist_notgermane" ~ "Bad for District, Not Germane",
                             rowname == "T_baddealdist_germane" ~ "Bad for District, Germane",
                             rowname == "T_notfairshare_notgermane" ~ "Fair Share, Not Germane",
                             rowname == "T_notfairshare_germane" ~ "Fair Share, Germane"))

make_kable(tablea22, col.names = c("Variable", "Mean", "SD", "Min", "Max"),
           caption = "Summary Statistics for Experiment 2")




# Table D2 ---------------------------------------------------------------------

lma23proj_dist <- lm_robust(project_thermometer ~  distabovezero + distworsecity +
                              city_pc + dist_pc + fielding + T_generic +
                              T_baddealdist_notgermane  +
                              T_notfairshare_notgermane + T_notfairshare_germane, data = df2,
                            subset = dist_pc >= 0)
lma23approve_dist <- lm_robust(approve_num ~  distabovezero + distworsecity +
                                 city_pc + dist_pc + fielding + T_generic +
                                 T_baddealdist_notgermane  +
                                 T_notfairshare_notgermane + T_notfairshare_germane, data = df2,
                               subset = dist_pc >= 0)
lma23incum_dist <- lm_robust(councilor_thermometer ~ distabovezero + distworsecity +
                               city_pc + dist_pc + fielding + T_generic +
                               T_baddealdist_notgermane  +
                               T_notfairshare_notgermane + T_notfairshare_germane, data = df2,
                             subset = dist_pc >= 0)
lma23vote_dist <- lm_robust(vote_num ~ distabovezero + distworsecity +
                              city_pc + dist_pc + fielding + T_generic +
                              T_baddealdist_notgermane  +
                              T_notfairshare_notgermane + T_notfairshare_germane, data = df2,
                            subset = dist_pc >= 0)
make_texreg(name = "tablea23",
            model = list(lma23proj_dist, lma23approve_dist, lma23incum_dist, lma23vote_dist),
            custom.model.names = c("\\makecell{Project\\\\Evaluation\\\\(0 to 100)}", 
                                   "\\makecell{Approval of\\\\Project\\\\(1 to 5)}", 
                                   "\\makecell{Incumbent\\\\Evaluation\\\\(0 to 100)}", 
                                   "\\makecell{Vote for Incumbent\\\\vs. Challenger\\\\(1 to 5)}"),
            custom.coef.map = list("distabovezero" = "District Benefits (District > 0)",
                                   "dist_pc" = "District Returns Per Capita",
                                   "distworsecity" = "District Worse Off than City",
                                   "city_pc" = "City Returns Per Capita",
                                   "fielding" = "Second Sample",
                                   "T_generic" = "Generic Critique",
                                   "T_baddealdist_notgermane" = "Bad Deal Critique (Not Germane)",
                                   "T_notfairshare_notgermane" = "Fair Share Critique (Not Germane)",
                                   "T_notfairshare_germane" = "Fair Share Critique (Germane)",
                                   "(Intercept)" = "Constant"),
            custom.note = paste("\\item[\\hspace{-5mm}] %stars.",
                                "\\item[\\hspace{-5mm}] \\textit{Note:} Dependent variables are listed in each column. Models estimated using ordinary least squares regression, 
       with standard errors clustered by respondent."),
            caption = "Effect of Challenger Criticisms on Evaluations when District At Least Breaks Even, Experiment 2",
            resize_width = TRUE, multiple.tasks = "normal")






# (Kable) Table D3 -------------------------------------------------------

tablea24 <- full_lrtest(df2)
make_kable(tablea24, col.names = c("", "$\\chi^2$",
                                   "$p$-value"),
           caption = "Likelihood Ratio Tests, Experiment 2",
           footnote = "\\textit{Note:} We conduct a multinomial logistic regression of the treatment variable
               on covariates party identification, household income, race, education,
               and region of the survey taker to test if our treatments of city and 
               district net return were associated with any key covariates. We then conducted
               a likelihood ratio test of this model with the null model for each treatment variable outcome.
               Each column reflects the likelihood ratio difference between a saturated and null model,
               as well as the $p$-value for rejecting the null hypothesis. Since we fail to reject
               the null for both variables, we are unable to say there is a significant
               association between our treatment variables and these covariates.")


# *Figures ----------------------------------------------------------------

# Figure E1 ---------------------------------------------------------------
lma1 <- lm_robust(eval_incumbent ~ paydist_cat + paycity_cat + task, 
                  data = df1, cluster = ResponseId)

KMG <- function(x){
  ifelse(x == 0, 0,
         ifelse(abs(x) < 1000000, paste0(x/1000, "k"),
                paste0(x/1000000, "m")))
}


dffa1 <- tidy(lma1) %>% 
  select(term, estimate, conf.low, conf.high) %>% 
  mutate(facet = case_when(str_detect(term, "dist") ~ "District Net Returns",
                           str_detect(term, "city") ~ "City Net Returns",
                           TRUE ~ NA_character_)) %>% 
  filter(!is.na(facet)) %>% 
  mutate(term = parse_number(term)) %>% 
  left_join(rbind(df1 %>% 
                    select(city, paycity_cat) %>%
                    rename(term = paycity_cat,
                           label = city) %>% 
                    mutate(facet = "City Net Returns",
                           term = as.numeric(term)), 
                  df1 %>% 
                    select(dist, paydist_cat) %>% 
                    rename(term = paydist_cat,
                           label = dist) %>% 
                    mutate(facet = "District Net Returns",
                           term = as.numeric(term))) %>% 
              distinct(label, facet, .keep_all = TRUE)) %>% 
  mutate(label = KMG(label),
         label_city = case_when((facet == "City Net Returns") ~ label),
         label_dist = case_when((facet == "District Net Returns") ~ label),
         label_city = fct_inorder(label_city),
         label_dist = fct_inorder(label_dist))


for(i in 1:2){
  type <- c('District Net Returns', 'City Net Returns')
  str_to_lower(str_sub(type[2], 1, 4))
  
  assign(paste0("p", i),
         ggplot(subset(dffa1, facet == 
                         c('District Net Returns', 'City Net Returns')[i])) +
           geom_pointrange(aes_string(x = paste0("label_", str_to_lower(str_sub(type[i], 1, 4))),
                                      y = "estimate", ymin = "conf.low", ymax = "conf.high")) +
           scale_y_continuous(breaks = c(-0.2, 0, 0.2, 0.4, 0.6, 0.8, 1.0),
                              limits = c(-.2, 1)) +
           facet_wrap(~facet) +
           geom_hline(yintercept = 0, linetype = "dotted") + 
           labs(x = "Net Returns in $", y = "Coefficient") + 
           theme_few())
  if(i == 2){
    p2 <- p2 + theme_few() + theme(axis.text.x = element_text(size = 6))
  }
  
}

gfa1 <- grid.arrange(p1,p2)

ggsave(paste0(save_path, "figurea1.png"),
       plot = gfa1, width = 6, height = 8)

# Figure E2 ---------------------------------------------------------------

lmfa2 <- lm_robust(eval_incumbent ~ paydist_cat + paycity_cat + task, 
                   data = df1R1, cluster = ResponseId)

dffa2 <- tidy(lmfa2) %>% 
  select(term, estimate, conf.low, conf.high) %>% 
  mutate(facet = case_when(str_detect(term, "dist") ~ "District Net Returns",
                           str_detect(term, "city") ~ "City Net Returns",
                           TRUE ~ NA_character_)) %>% 
  filter(!is.na(facet)) %>% 
  mutate(term = parse_number(term)) %>% 
  left_join(rbind(df1R1 %>% 
                    select(city, paycity_cat) %>%
                    rename(term = paycity_cat,
                           label = city) %>% 
                    mutate(facet = "City Net Returns",
                           term = as.numeric(term)), 
                  df1R1 %>% 
                    select(dist, paydist_cat) %>% 
                    rename(term = paydist_cat,
                           label = dist) %>% 
                    mutate(facet = "District Net Returns",
                           term = as.numeric(term))) %>% 
              distinct(label, facet, .keep_all = TRUE)) %>% 
  mutate(label = KMG(label),
         label_city = case_when((facet == "City Net Returns") ~ label),
         label_dist = case_when((facet == "District Net Returns") ~ label),
         label_city = fct_inorder(label_city),
         label_dist = fct_inorder(label_dist))

for(i in 1:2){
  type <- c('District Net Returns', 'City Net Returns')
  str_to_lower(str_sub(type[2], 1, 4))
  
  assign(paste0("p", i),
         ggplot(subset(dffa2, facet == 
                         c('District Net Returns', 'City Net Returns')[i])) +
           geom_pointrange(aes_string(x = paste0("label_", str_to_lower(str_sub(type[i], 1, 4))),
                                      y = "estimate", ymin = "conf.low", ymax = "conf.high")) +
           scale_y_continuous(breaks = c(-0.2, 0, 0.2, 0.4, 0.6, 0.8, 1.0),
                              limits = c(-.2, 1)) +
           facet_wrap(~facet) +
           geom_hline(yintercept = 0, linetype = "dotted") + 
           labs(x = "Net Returns in $", y = "Coefficient") + 
           theme_few())
  if(i == 2){
    p2 <- p2 + theme_few() + theme(axis.text.x = element_text(size = 6))
  }
  
}

gfa2 <- grid.arrange(p1,p2)

ggsave(paste0(save_path, "figurea2.png"),
       plot = gfa2, width = 6, height = 8)

# Figure E3 ----------------------------------------------------------------

xlabels_cat3dist <- c("-$50", "-$40", "-$30", "-$25", "-$20", "-$15", "-$10",
                      "-$5", "-$1", "-$0.1", "$0", "$0.1", "$1", "$5", "$10",
                      "$15", "$20", "$25", "$30", "$40", "$50")

breaks = c(-50, -40, -30, -25, -20, -15, -10, -5, -1,
           -0.1, 0, 0.1, 1, 5, 10, 15, 20, 25, 30, 40, 50)

# Step 2 - Visualize New Data Frame
df_cat3dist_count <- df1 %>% 
  mutate(cat3dist = case_when(dist > 0 ~ "Net Gain",
                              dist == 0 ~ "Break Even",
                              dist < 0 ~ "Net Loss")) %>% 
  group_by(cat3dist, city_pc) %>% 
  count()

df_cat3dist <- df1 %>% 
  mutate(cat3dist = case_when(dist > 0 ~ "Net Gain",
                              dist == 0 ~ "Break Even",
                              dist < 0 ~ "Net Loss")) %>% 
  group_by(cat3dist, city_pc) %>% 
  summarize(meaneval = mean(eval_incumbent, na.rm = TRUE)) %>% 
  left_join(df_cat3dist_count) %>% 
  mutate(n = n^1.5,
         cat3dist = factor(cat3dist,
                           levels = c("Net Gain", "Break Even", "Net Loss")),
         city_pc = factor(city_pc))%>% 
  rename(sizing = n,
         `District Return` = cat3dist)

options(ggplot2.discrete.colour= c("blue", "black", "red"))


gfa3 <- ggplot(df_cat3dist, aes(x = city_pc, y = meaneval, col = `District Return`)) +
  geom_line(aes(group = `District Return`)) +
  geom_point(aes(size = sizing, shape = `District Return`)) +
  scale_shape_manual(values = c("\u002b","\u0030",  "\u002d")) + #change values to plus, zero, minus
  ylim(-1,1) + 
  scale_size(range = c(0, 5), guide = 'none') +
  scale_x_discrete(breaks = breaks,
                   labels = xlabels_cat3dist) + 
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(x = "City Per Capita Return (Non-Linear Scale)", 
       y = "Evaluation of Incumbent") +
  theme_few() +
  theme(
    legend.position = c(0.9, 0.12),
    plot.title = element_text(hjust = 0.5)) +
  annotate("rect", xmin=10, xmax=11, ymin=-Inf, ymax=Inf, alpha=0.5, fill="gray") +
  geom_fit_text(min.size = 2.5,
                data = data.frame(xmin = 10, xmax = 11,
                                  label = "Expected\nDiscontinuity\nInterval"),
                aes(xmin = xmin, xmax = xmax, ymin = 0.5, ymax = 0.8, label = label,
                    fontface = "bold"),
                inherit.aes = FALSE)
gfa3


ggsave(paste0(save_path, "3categorydist.png"),
       plot = gfa3, width = 8, height = 6)
# Figure E4 ----------------------------------------------------------------


xlabels_cat3city <- c("-$50", "-$25", "-$10", "-$5", "-$1",
                      "$0", "$1", "$5", "$10", "$25", "$50")

breaks = c("-1250000","-625000", "-250000", "-125000", 
           "-25000", "0", "25000", "125000", "250000",
           "625000", "1250000")

# Step 2 - Visualize New Data Frame
df1R1_cat3city_count <- df1R1 %>% 
  mutate(cat3city = case_when(city > 0 ~ "Net Gain",
                              city == 0 ~ "Break Even",
                              city < 0 ~ "Net Loss")) %>% 
  group_by(cat3city, dist_pc) %>% 
  count()

df1R1_cat3city <- df1R1 %>% 
  mutate(cat3city = case_when(city > 0 ~ "Net Gain",
                              city == 0 ~ "Break Even",
                              city < 0 ~ "Net Loss")) %>% 
  group_by(cat3city, dist_pc) %>% 
  summarize(meaneval = mean(eval_incumbent, na.rm = TRUE)) %>% 
  left_join(df1R1_cat3city_count) %>% 
  mutate(n = n^2,
         cat3city = factor(cat3city,
                           levels = c("Net Gain", "Break Even", "Net Loss")),
         dist_pc = factor(dist_pc))%>% 
  rename(sizing = n,
         `City Return` = cat3city)

options(ggplot2.discrete.colour= c("blue", "black", "red"))


gfa4 <- ggplot(df1R1_cat3city, aes(x = dist_pc, y = meaneval, col = `City Return`)) +
  geom_line(aes(group = `City Return`)) +
  geom_point(aes(size = sizing, shape = `City Return`)) +
  scale_shape_manual(values = c("\u002b","\u0030",  "\u002d")) + #change values to plus, zero, minus
  ylim(-1,1) + 
  scale_size(range = c(0, 20), guide = 'none') +
  scale_x_discrete(breaks = c(-50, -25, -10, -5, -1, 0, 1, 5, 10, 25, 50),
                   labels = xlabels_cat3city) + 
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(x = "District Per Capita Return (Non-Linear Scale)", 
       y = "Evaluation of Incumbent") +
  theme_few() +
  theme(
    legend.position = c(0.9, 0.12),
    plot.title = element_text(hjust = 0.5)) +
  annotate("rect", xmin=5, xmax=6, ymin=-Inf, ymax=Inf, alpha=0.5, fill="gray") +
  geom_fit_text(
    data = data.frame(xmin = 5, xmax = 6, 
                      label = "Expected\nDiscontinuity\nInterval"),
    aes(xmin = xmin, xmax = xmax, ymin = 0.5, ymax = 0.8, label = label,
        fontface = "bold"),
    inherit.aes = FALSE)
gfa4


ggsave(paste0(save_path, "figurea4.png"),
       plot = gfa4, width = 8, height = 6)



# Figure E5 ---------------------------------------------------------------


xlabels_cat3city <- c("-$55", "-$50", "-$45", "-$30", "-$25", "-$20", "-$15", "-$10", "-$5", "-$2", "-$0.5", "-$0.1",
                      "$0", "$0.1", "$0.5", "$2", "$5", "$10", "$15", "$20", "$25", "$30", "$45", "$50",  "$55")

# Step 2 - Visualize New Data Frame
df_replication2_count <- df1R2 %>% 
  group_by(aggregate_order, dist_pc) %>% 
  count()

df_replication2 <- df1R2 %>% 
  group_by(aggregate_order, dist_pc) %>% 
  summarize(meaneval = mean(eval_incumbent, na.rm = TRUE)) %>% 
  left_join(df_replication2_count) %>% 
  mutate(n = n^2,
         dist_pc = factor(dist_pc)) %>% 
  rename(sizing = n)

options(ggplot2.discrete.colour= c("#F8766D", "#00BFC4"))

ggplot(df_replication2, aes(x = dist_pc, y = meaneval, col = aggregate_order)) +
  geom_line(aes(group = aggregate_order, linetype = aggregate_order)) +
  geom_point(aes(size = sizing))+
  #scale_shape_manual(values = c("\u002b","\u0030",  "\u002d")) + #change values to plus, zero, minus
  ylim(-1,1) + 
  scale_linetype_manual(values=c("solid", "twodash")
                          )+
  scale_size(range = c(0, 10), guide = 'none')+
  scale_x_discrete(breaks = c(-55, -50, -45, -30, -25, -20, -15, -10, -5, -2, -.5, -0.1, 0, 0.1, 0.5, 2, 5, 10, 15, 20, 25, 30, 45, 50, 55),
                   labels = xlabels_cat3city)+
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(x = "District Per Capita Return (Non-Linear Scale)", 
       y = "Evaluation of Incumbent",
       linetype = "Aggregate\nJurisdiction",
       col = "Aggregate\nJurisdiction") +
  theme_few()+
  theme(
    legend.position = c(0.9, 0.14),
    plot.title = element_text(hjust = 0.5)) +
  annotate("rect", xmin=12, xmax=13, ymin=-Inf, ymax=Inf, alpha=0.5, fill="gray") +
  geom_fit_text(min.size = 2.5,
                data = data.frame(xmin = 12, xmax = 13,
                                  label = "Expected\nDiscontinuity\nInterval"),
                aes(xmin = xmin, xmax = xmax, ymin = 0.5, ymax = 0.8, label = label,
                    fontface = "bold"),
                inherit.aes = FALSE)

ggsave(paste0(save_path, "figure_e5.png"), width = 9, height = 6)



